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Abstract. The description of muon propagation through 
large depths of matter, based on a concept of the correction 
factor, is proposed. The results of Monte-Carlo calculations 
of this correction factor are presented. The parametrizations 
for continuous energy loss coefficients, valid in the broad 
interval of muon energies, and for the correction factor are 
given. The concrete calculations for pure water are pre- 
sented. 



1 Introduction 

The existing methods of calculations of muon propagation 
through thick layers of matter (basically for a standard rock) 
are largely presented in the literature. There are two com- 
pletely different approaches in such calculations. In first 
(historically) works on the subject of muon propagation the 
integro-differential kinetic equation for muon flux had been 
formulated and approximately solved using semianalytical 
methods. We mention here one of the pioneering works of 
this approach: (Zatsepin and Mikhalchi, 1962), and one of 
the most recent papers: (Naumov, Sinegovsky, and Bugaev, 
1993) (this latter paper contains also the large bibliography). 
Works of the second approach use Monte Carlo (MC) tech- 
nique for propagation studies: the main element in such 
works is simulation of sequences of free flights and inter- 
actions during muon passage through a medium. One of the 
most early works using this method is (Hayman, Palmer, and 
Wolfendale, 1963), the most recent one is the work of the 
present authors (Sokalski, Bugaev, and Klimushin, 2000), 
containing also many bibliographic references. 

2 Method 

The integral flux of muons in matter with energies above cut- 
off Ef expected at slant depth i? at a zenith angle of 9 taking 
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into account fluctuating character of muon losses is conven- 
tionally described by 

/•OO 

Ifi{>Ef,R,0)^ / P{Eo,R,>Ef)D{Eo,e)dEo, (1) 
Jo 

where P{Eo, R,> Ef) is the probability that a muon, hav- 
ing starting energy Eq, after passing of path R will sur- 
vive with final energy above cut-off Ef and D{Eo,9) = 
dN{Eo, 9)/dEo is a sea level differential angular spectrum. 
In a case of a flat surface the slant depth is expressed by 
R = h/ cos 6, where his a vertical depth below the surface. 

With the assumption of continuous energy loss rate of 
muon in matter, L{E) = -dE/dx = a{E) + h{E)E, the 
integral flux Id is derived from the equation 

/•OO 

Ici{>Ef,R,e)= D{Eo,0)dEo, (2) 

where the low Umit of the integral. Eg, is the value of starting 
energy Eq which results, after passing of path R, in the final 
energy Ef. This value is derived from the solution of the 
integral equation 

/■^^ dE 

The method proposed by us for a calculation of real in- 
tegral muon flux allowing for loss fluctuations consists in 
the following. As is shown in (Klimushin, Bugaev, and 
Sokalski, 2000), if we approximate the muon energy loss 
function L{E) by a linear energy dependence, 

L{E) =a + l3E (4) 

(where a and (3 are independent on energy), and if we 
parametrize the sea level muon spectrum by dependencies 
of the type 
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Fig. 1. The b-terms of muon radiative energy losses in pure water 
(as MUM code output). Losses due to e+e" pair production bp 
(curve (a)), bremsstrahlung bt (b), and photonuclear interaction 6„ 
(c) are shown. Curve (d) is the sum of all b-terms. 



Fig. 2. Muon energy losses in pure water as a function of energy 
E (as MUM code output). Curve (a) is loss due to ionization a{E) 
and (b) is total loss due to all radiative processes b{E)E. Curve (c) 
describes total energy losses a{E) + b{E)E. 



then the integral flux Id (Eq.6) can be expressed analytically. 
Further, to calculate a real Ifi we must know the correction 
factor defined by the ratio 



Cf{>Ef,R,e) 



i,i{>Ef,R,e) 
ifii>Ef,R,ey 



(5) 



In principle, this factor can be calculated using known codes 
for muon propagation through matter. In this work we apply 
for this aim the MUM code described in previous work of 
the present authors (Sokalski, Bugaev, and Klimushin, 2000). 
As we will see, this factor, as a function of muon energy 
and depth, can be easily parametrized by simple expressions. 
As a result, the real flux //; is conveniently expressed by 
the formula, containing continuous energy loss coefficients 
a and (3. 



3 Continuous energy losses 

For the description of continuous energy losses of muon 
in water needed for solving the integral equation (3) we 
have done parametrizations based on output data from MUM 
code. All formulas for cross sections used in MUM code for 
muon energy loss computations are described in details in 
Ref. (Sokalski, Bugaev, and Klimushin, 2000). 

The energy loss rate per unit of path x is given, conven- 
tionally, by the expression L{E) = —dE/dx = a(E) + 
b{E)E, where a{E) is loss due to ionization and b{E) = 
bp{E) + bb{E) + bn{E) is the sum of coefficients for aU ra- 
diative processes: e+e" pair production (&p), bremsstrahlung 
(5b) and photonuclear interaction (6„). 



Ionization loss was taken as composition of 2 processes 
a{E) = ac{E) + ae{E), where Uc is classic ionization cal- 
culated using Bethe-Bloch formula and Ue results from e- 
diagrams for bremsstrahlung being treated as a part of ion- 
ization process following (Kelner, Kokoulun, and Petrukhin, 
1997) (7-quantum is emitted by atomic electron). Taking 
into account the last process leads to 1.8 %, 3.4 %, 5.5 % in- 
crease of ionization loss for 100 GeV, 1 TeV, 10 TeV muons, 
correspondingly. The approximation formula for Cc is given 
by 



ac{E) = flco + Oci In 



W = 



w„ 



E 



(6) 



l + m2/(2me£)' 

where Wmax is maximum energy transferable to the electron 
and m^, me are the rest masses of muon and electron. The 
set of cofficients, in units of (10~^GeVcm^g~^), 

flco = 2.106, ttci = 0.0950 for E<45 GeV, 
= 2.163, ttci = 0.0853 for E > 45 GeV, 

gives the error for parametrization (6) smaller than 0.2 % for 
(1-10^) GeV range. For Oe the following polynomial ap- 
proximation {E in units of (GeV)): 

aeiE) =3.54 3.785 ln£-|-1.151n^i; 
-1-0.0615 In^i; (10-^GeVcm^g-^), (7) 

has the error < 0.2 % for (50-10^) GeV range. Finally, 
the sum of (6) and (7) has the error of a{E) approximation 
smaller than 0.2 % for (1-10*) GeV range. 
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Table 1. Coefficients bij (10 ®cm^g ^) of the fitting formula (8) for b-terms of muon energy losses in water. Maximum absolute values of 
relative errors are shown. 
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Radiative energy losses of muon in water has been cal- 
culated using the cross sections from the following works: 
(Andreev, Bezrukov, and Bugaev, 1994) for bremsstrahlung, 
(Kelner, Kokoulin, and Petrukhin, 1999) for e+e^ pair pro- 
duction, and (Bezrukov and Bugaev, 1980) for photonuclear 
interaction. Fig. 1 sketches corresponding b-terms of radia- 
tive energy losses. 

In principle MUM code gives values for all b-terms up 
to 10^ GeV. Note that the logarithmic rise with an energy 
of the photoabsorption cross section, used in the model of 
Ref. (Bezrukov and Bugaev, 1980) and resulting in corre- 
sponding increase of the 6„-term, is experimentally proved 
already up to the photon energy ^ 10'' GeV. 

For parametrization of b-terms we have used the same 
functional form as in Ref. (Bezrukov and Bugaev, 1981) but 
increased the power of polynomial to improve accuracy and 
enlarge the range of application up to 100 PeV without divi- 
sion on energy subintervals {E in units of (GeV)): 
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bi{E) = ^ hij In-' E, where i = p, b, n. (8) 

Corresponding coefficients of this decomposition are col- 
lected in Tab. 1. This fit works with typical errors ±(0.2- 
0.4)% within (50-7 X 10^) GeV for all b-terms. 
The sum b{E) of fits for b-terms is valid for the energy range 
(50-10^) GeV with the relative error < 0.5 %. 

Fig. 2 shows muon energy loss rate in water. Note that the 
radiative losses account for approximately 5 % of the value 
of total losses at muon energy 70 GeV, 50 % at 1 TeV and 
dominate as 95 % at 20 TeV. The sum of all parametriza- 
tions (6), (7) and (8) results in the description of total energy 
losses L{E) with excellent accuracy < 0.3 % with varying 
sign of error for energy range (1-10*) GeV mainly because 
of b-term error compensations (even for region (1-50) GeV). 
Thanks to this fact just this parametrization of L{E) has been 
used for the numerical computation of integral flux Id de- 
fined by (2) and (3). 

Further, we realized that the total losses L{E) may be de- 
scribed with an accuracy better than 2.5 % (with varying sign 
of error) for energy range (3-3 x 10^) GeV by 3-slope hnear 
fit 

L{E) =a + (3E, (9) 
where 

a = ao = 2.30, f3 = Po = 15.50 for E < Eqi = 30.0GeV, 



a = ai = 2.67,/3 = /3i = 3.40 

fori^oi <E <Ei2^ 35.3TeV, 
a = a2 = -6.50, /? = /32 = 3.66 for E > E12 = 35.3TeV, 

and a are in (IG^'^GeVcm^g^^) and (3 in (lO^'^cm^g^^). 
The energy losses expressed by this parametrizations have 
sense of effective ones, for example in energy region 30 
GeV-35 TeV the values a = 2.67 and /? = 3.40 represent 
effective energy losses due to ionization and radiative pro- 
cesses, correspondingly. The availability of the linear depen- 
dence (9) leads to the possibility to derive an analytical for- 
mula for underwater integral flux (see (Klimushin, Bugaev, 
andSokalski, 2001)). 



4 Correction factor 

The influence of fluctuations of muon losses in matter (mainly 
due to radiative processes) results in that the real integral flux 
//; is generally greater than Id calculated in the approxima- 
tion of continuous losses. In this work we propose to al- 
low for the influence of energy loss fluctuations on the value 
of angular flux in matter by means of correction factor ex- 
pressed by the ratio (5). In the assumption when the same 
differential cross sections are used for computing both the 
numerator and denominator of ratio (5), one may expect that 
this factor depends only weakly on sea level angular spec- 
trum. For numerical calculations of the correction factor we 
used total continuous energy losses defined by sum of (6), (7) 
and (8), and survival probabilities calculated by using MUM 
code (Klimushin, Bugaev, and Sokalski, 2000). The values 
of correction factors calculated for the same slant depth R at 
vertical direction and at zenith angle 9 differ weakly. It is 
iUustrated in Fig. 3, where one can see that C/(> Ef, R, 0°) 
differs from C/(> Ef, R, arccos /i/i?) maximum on 3.3 % 
for Ef >10 GeV at vertical depth /i of 1.15 km. It appears 
that with acceptable accuracy the correction factor depends 
on slant depth R only, rather than on R and 9 separately. 

The dependencies of correction factor on Ef and R, cal- 
culated for sea level spectrum given by expression (5) of 
Ref. (Klimushin, Bugaev, and Sokalski, 2001), represent the 
set of rather smooth curves (shown in Fig. 3) and it is possi- 
ble to approximate this factor by formula 

4 4 

Cf{> Ef, R, ^) = Cij logio Ef)R\ (10) 

i=0 j=0 
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Table 2. Coefficients dj of the fitting formula (10) for correction factor calculated for vertical sea level spectrum given by expression (5) of 
Ref. (Klimushin, Bugaev, and Sokalski, 2001). 



subscripts do cn ca ca cn 






6.3045 X 10" 


1 


6.6658 X 10" 


1 


-4.5138 X 10" 


1 


1.2441 X 10" 


i 


-1.1904 X 10" 


-2 


1 


2.0152 X 10" 


1 


-4.2990 X 10" 


1 


3.2532 X 10" 


1 


-1.0265 X 10" 


1 


1.0751 X 10" 


-2 


2 


-3.3419 X 10" 


2 


5.1833 X 10" 


2 


-3.9229 X 10" 


2 


1.2360 X 10" 


2 


-1.2911 X 10" 


-3 


3 


1.6365 X 10" 


3 


-2.3645 X 10" 


3 


1.7775 X 10" 


3 


-5.5495 X 10" 


4 


5.7557 X 10" 


-5 


4 


-2.6630 X 10" 


5 


3.7770 X 10" 


5 


-2.8207 X 10" 


5 


8.7275 X 10" 


6 


-8.9919 X 10" 


-7 




R (km) 

Fig. 3. Correction factor C/ as a function of slant depth R in pure 
water. The results obtained using sea level spectrum defined by ex- 
pression (5) of Ref. (Klimushin, Bugaev, and Sokalski, 2001) are 
given. Solid curves correspond to numerical calculations for ver- 
tical case 9 — 0° . Dashed curves describe the correction factor 
computed at vertical depth h of 1.15 km for various zenith angles 
as a function of slant depth defined by R = h/ cos 9. Both solid 
and dashed curves are shown for four values of cut-off energy Ef-. 
10 GeV, 100 GeV, 1 TeV, and 10 TeV, from top to bottom. 

Here cut-off energy Ef is expressed in (GeV) and slant depth 
R is in (km) with the coefficients collected in Table 2. 
When using (10) for cut-off energies Ef <10 GeV one 
should substitute value of Ef =10 GeV. 

Formula (10) can be applied for any geometrical shape of 
the surface. Right hand side of (10) depends on 9 because, 
generally, R = R{6). So, in the particular case of a flat sur- 
face the angular dependence of the correction factor appears, 
in our approximation, only through the relation R = h/ cos 6 
(where his a vertical depth). 

The accuracy of formula (10) for _Ej^=(l-100) GeV is bet- 
ter than ±2 % for slant depths R as large as 22 km and is 
not worse than ±3 % for Ef=l TeV up to i?=17 km and for 



Ef=lO TeV up to R=15 km. Fig. 3 shows that for Ef < 100 
GeV the total energy loss may be treated as quasi-continuous 
(at level of C/ > 0.9) only for slant depths R < 2.5 km but 
for Ef=lO TeV the fluctuations should be taken into account 
at level of 15 % already for slant depth as small as R=l km. 
For slant depths larger than 10 km Ifi and Id differ more 
than on a factor of 2. 

5 Conclusions 

We showed in this paper that the most complicate part of 
any calculation of muon integral flux at great depths (which 
is taking into account the energy loss fluctuations) can be 
reduced to a computation of the correction factor which de- 
pends only weakly on the sea level spectrum and is a smooth 
function of its arguments (muon cut off energy and a slant 
depth). The method described here permits to carry out 
fast analytical calculations of underwater and underground 
fluxes, as is shown in Ref. (Klimushin, Bugaev, and Sokalski, 
2001). These analytical calculations are greatly facilitated 
by the parametrizations of all continuous energy losses (us- 
ing the most recent expressions for muon interaction cross 
sections) proposed in the present paper. 
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